arXiv:1507.06957vl [cond-mat.soft] 22Jul2015 


Addendum to “Vector-based model of elastic bonds for simulation of granular solids 

Vitaly A. Kuzkin* 

Institute for Problems in Mechanical Engineering RAS, 

Saint Petersburg Polytechnical University 
(Dated: July 27, 2015) 

In our previous work [Phys. Rev. E 86, 05130 (2012)] the model for simulation of granular 
solids composed of bonded particles was presented. In the present Brief Report the method for 
validation of numerical implementation of the model is proposed. Four benchmark problems are 
solved analytically. Expressions for forces and torques acting between two bonded particles in the 
case of tension, shear, bending, and torsion are presented. Relations between parameters of the 
model, bond stiffnesses, and mechanical properties of bonding material are derived. Additionally, 
the expression for potential energy of the bond is substantially simplified. It is shown that, in spite 
of simplicity, the model is applicable to bonds with arbitrary length/thickness ratio. 
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I. INTRODUCTION 

In our previous work [1] the model for simulation of 
solids composed of bonded rigid-body particles was de¬ 
rived. The bonds either represent additional glue-like 
material connecting particles [2-5] or appear as a result 
of coarse-graining of macromolecules [6], nanotube-based 
materials [7], etc. In both cases the bond causes forces 
and torques acting on the bonded particles. In paper [1] 
the potential energy of the bond is represented as a func¬ 
tion of vectors, rigidly connected with the particles. Cor¬ 
responding expressions for forces and torques are derived. 
Several examples of simulation of rod-like and shell-like 
granular solids are presented. However in the given ex¬ 
amples the behavior of the system is quite complicated. 
As a result, validation of numerical implementation of 
the bond model is not straightforward [8]. 

In the present Brief Report simple approach for valida¬ 
tion is presented. Four test problems for a system of two 
bonded particles are solved analytically. The resulting 
expressions for forces and torques can be used for valida¬ 
tion of computer codes. Additionally, the expression for 
potential energy of the bond is substantially simplified. 
Relations between parameters of the model, characteris¬ 
tics of bonding material, and stiffnesses of the bond are 
derived. In spite of its simplicity, the model has the same 
advantage as the original model derived in paper [1]: it 
allows to fit any values of the bond stiffnesses exactly. 


II. SIMPLIFIED MODEL OF A SINGLE BOND 

Consider simple model of an elastic bond in a granular 
solid composed of particles. In general, every particle can 
be bonded with any number of neighbors. The behavior 
of the bonds is assumed to be independent, i.e. pairwise 
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interactions are considered. Therefore for simplicity two 
bonded particles i and j are considered below. The ex¬ 
pressions for forces and torques acting on the particles i 
and j are derived. 

Assume that the bond connects two points that belong 
to the particles. The points lie on the line connecting the 
particles’ centers in the initial (undeformed) state. For 
example, the points can coincide with particles centers. 
Let us denote distances from the points to particles’ cen¬ 
ters of mass as Ri, Rj respectively. For example, in the 
case shown in figure 1, the points lie on particles’ surfaces 
and values Ri, Rj coincide with particles’ radii. 

The potential energy of the bond is represented 
as a function of orthogonal unit vectors n,;i,rq 2 ,rq 3 
and nji,rij 2 ,nj 3 , rigidly connected with particles i and 
j respectively. The first index corresponds to particle 
number, the second index corresponds to vector’s num¬ 
ber (see figure 1). In the undeformed state the following 



FIG. 1: Two bonded particles in the undeformed state (left) 
and deformed state (right). Here and below a is an equilib¬ 
rium distance. 

relations are satisfied: 

fly 1 — Gij , Ui2 4R3 Uy3 . (1) 

where = r^/r^, r^ = rj r t . In paper [1] it is 
shown that the potential energy of the bond is a func¬ 
tion of vector D, y = f + R 3 nj-| — RpM and vec¬ 
tors nik,nj m ,k,m = 1,2,3. Vector D,j connects the 
points defined by vectors r i + Rinn, rj + RjUji (see fig¬ 
ure 1). In the case Ri = Rj = 0 the bond connects par¬ 
ticle centers. In the present Brief Report the following 
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expression for potential energy of the bond is proposed: 
B B 

U = y (Dij - a) 2 + y (nji - n a ) • dij + B 3 n a • n jl 
-y'(n i2 ■ n j2 + n i3 • n j3 ). 

(2) 

where Djj — D jj . d / i; — F) jj/ij ■ B\ . B 2 , B ' 3 , -64 are 
parameters of the model. In the following section it is 
shown that parameters B are related to stiffnesses of 
the bond. Note that the expression (2) is significantly 
simpler than the one proposed in paper [1], The main 
difference is in the last term. The last term in formula (2) 
contributes to both bending and torsion of the bond. In 
paper [ 1 ] the potential energy is designed in such a way 
that at small deformations bending and torsion are de¬ 
scribed by two independent terms. This “decomposition” 
leads to unnecessarily complicated expression for poten¬ 
tial energy. 

Forces F,,- = — Fj* and torques M,,-. M,;,- acting be¬ 
tween the particles have the form: 

B 2 

F ij = — a)dij + 2 j) ( n .A — n *i) ' — ;)• 

My = Riia.il x F ij - ij x na + M TB , 

Mji = Rj nj 1 x F ji + yd ij x n,-| - M TB , 

B 4 

M TB = B 3 TOji x n.;-| - — (n ,2 x n l2 + n j3 x n i3 ), 

(3) 

where E is a unit tensor. If the bond connects particle 
centers, the expressions (3) take form: 

B 2 

F ij B[ (j'ij (l) r , j • — (11,1 II,1 j * (E GijGij ), 

Ifij 

My = -y eij x n,;! + M TB , 
x n ji - M TB , 

(4) 

where = r ij/rij, M TB is defined by equation (3). 

Note that the forces and torques (3), (4) are defined by 
instantaneous positions and orientations of the particles. 


III. BENCHMARKS: STRETCHING, SHEAR, 
BENDING, AND TORSION OF THE BOND 

The behavior of a solid composed of bonded particles 
interacting by forces and torques (3), in general, is quite 
complicated. Therefore validation of corresponding com¬ 
puter codes is not straightforward. In the present sec¬ 
tion, simple approach for validation is proposed. Forces 
and torques acting on two bonded particles are calculated 
in the case of stretching, shear, bending, and torsion of 
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tension shear bending torsion 

FIG. 2: Deformations of the bond and corresponding orien¬ 
tation of vectors, connected with the particles. 


the bond (see figure 2). The resulting expressions can 
be used as benchmarks for validation of computer im¬ 
plementation of the model. Additionally, the relations 
between parameters of the model and stiffnesses of the 
bond are derived. For simplicity it is assumed that the 
bond connects particle centers (/?,; = Rj = 0 ). 


A. Stretching and shear of the bond 


Consider pure stretching of the bond connecting parti¬ 
cles i and j. Stretching is carried out along the vector r ij. 
In this case forces and torques (3) have the form 

F ^ = Bi(rij - a)e zj , M ti = M^ = 0. (5) 


One can see that in the case of pure stretching the be¬ 
havior of the bond is identical to the behavior of a linear 
spring with stiffness B 3 . Therefore longitudinal stiffness 
of the bond ca is equal to B\. 

Consider pure shear in the two particle system. As¬ 
sume that position of particle i is fixed and particle j 
has displacement wk, where the unit vector k is orthog¬ 
onal to the initial direction of the bond eo- In this case 
forces and torques acting between the particles are the 
following: 


B 2 


F ij —B \ {j'ij o ) (eo cy • e,y 0,4 (, 

Vij ( 6 ) 

B 2 

Nfq* = Nfy i = 2 x ®o* 


Rewriting formulae ( 6 ) using the geometrical relations 


a u u 

e 0 * &ij — 5 &ij ’ k — , &ij x ©o — k x ©05 

f'ij T H T j/ 


' ij 
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Vij = 'Jo? + u 2 , 
yields 

F ^ • e 0 = Bia ( 1 - 
F ^ ■ k = Bin I 1 — 


(7) 


B 2 u 2 


Va 2 +u 2 J (a 2 +u 2 )2 
B 2 ua 


Jd 2 + u 2 J (a 2 +n 2 )2 


( 8 ) 


M i5 - = =- B ' 2 " 

3 3 2 Va 2 + u 2 


k x eg. 


Linearization of the second formula from ( 8 ) with respect 
to u yields Fjj-k = c^u, where cjj = B 2 /a 2 is the shear 
stiffness of the bond. 
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Thus parameters of the model B i and B -2 are propor¬ 
tional to longitudinal, ca, and shear, eg, stiffnesses of 
the bond respectively. 


B. Bending and torsion of the bond 


Consider pure bending of the bond. Assume that po¬ 
sitions of particles i and j are fixed. The particles are 
rotated in opposite directions by angle tp around vec¬ 
tor n *2 = n-j- 2 - In this case the forces vanish and the 
torques are equal with opposite sign: 


.\r„ - m , - • 



sinyj + B 3 sin2<p + 



n*2, 


A. Long bonds 


Mechanical behavior of relatively long bonds can be de¬ 
scribed by the beam theory [9]. Let us derive the relation 
between parameters of the model Bk and characteristics 
of massless Timoshenko beam connecting particles. As¬ 
sume that the beam has equilibrium length a, constant 
cross section, and isotropic bending stiffness. Longitu¬ 
dinal, shear, bending, and torsional stiffnesses of Timo¬ 
shenko beam are derived in paper [ 1 ]: 



ES 


CA = 

a 

EJ 

CD = 

c B = 

> 

a 

ct = 


12 kEJS 

a{nSa 2 + 24,7(1 + v)) ’ 
GJ p 
a 


( 12 ) 


F ij — 0. 

(9) 

In the case of small rotations of the particles = 
—Mjj« —2cB<^rij2, cb = B 2 /4 + B 3 + Bi/2, where c B is 
a bending stiffness of the bond. 

Consider torsion of the bond. Assume that position 
and orientation of the particle i is fixed and particle j is 
rotated around r^ by angle ip. In this case forces and 
torques (3) acting on the particles have the form 

B/± 

M ij = ~— {p.j2 x 11,2 +11,3 x n l3 ) = -B4smpeij, 
M ji = -M ij , Fij = 0. 

( 10 ) 

It is seen that torsional stiffness is equal to parameter B 4 
of the model, i.e. ct = B 4 . 

Thus stiffnesses of the bond are related to parameters 
of the model (2) by the following simple formulas: 


where E, G , v are Young’s modulus, shear modulus, and 
Poisson ratio of the bonding material; S, J, J v are cross 
section area, moment of inertia, and polar moment of 
inertia of the cross section respectively; k is a dimension¬ 
less shear coefficient [9]. Then formulas (11) and (12) 
yield the relation between parameters of the model and 
characteristics of Timoshenko beam: 

„ ES „ \2naEJS „ GJ V 

Bi = —, Bz ~ c 2 I o/i rn . B ± = —, 

a kSci z -)- 24j (1 -t- v) a 

EJ B 2 B 4 

B 3 = --. 

a 4 2 

(13) 

In the limit n —> oo Timoshenko beam is equivalent to 
Bernoulli-Euler beam. Corresponding relation between 
the parameters has the form: 



B 2 


\2EJ 

a 


2 EJ GJ p 
a 2 a 


B 2 B 2 ., B 4 

ca = B 1 , c D = — c B = —- + B 3 + Si—, c T = B 4 . 
a z 4 2 

(11 ) 

It is seen that by choosing parameters of the model any 
values of longitudinal, shear, bending, and torsional stiff¬ 
nesses of the bond can be fitted. Therefore the model is 
applicable to bonds with arbitrary length/thickness ra¬ 
tio. Thus in spite of the simplification, the model pre¬ 
sented above has the same advantage as the original one 
developed in paper [ 1 ]. 

The expressions (5), ( 8 ), (9), (10) can be used for val¬ 
idation of numerical implementation of the bond model 
proposed in the present paper. 


(14) 

If the parameters are determined by formula (14), then 
under small deformations the model is equivalent to a 
Bernoulli-Euler beam connecting particles. 

Thus formulas (13), (14) can be used for calibration of 
the model ( 2 ) in the case of relatively long bonds that 
can be approximated by Bernoulli-Euler or Timoshenko 
beam. 


B. Short bonds 


IV. RELATIONS BETWEEN PARAMETERS OF 
THE MODEL AND MECHANICAL 
CHARACTERISTICS OF THE BOND 

In the present section the relations between parame¬ 
ters of the model and mechanical properties of bonding 
material are derived. 


In materials composed of glued particles the bonds are 
usually short [5]. For example, the length/thickness ra¬ 
tio of the bonds in the material considered in paper [5] 
is much less than unity. In this case the beam models 
described above are inaccurate. Therefore in the present 
section an alternative model [ 1 ] is used for calibration of 
parameters B 3 , _B 2 , B 3 , B 4 . Longitudinal, shear, bend¬ 
ing, and torsional stiffnesses of a short bond are related 
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to characteristics of bonding material as follows [1]: 


(1 - v) ES 
(1 + i/)(l — 2v) a 


(1 — iz) EJ \_jup 

cb = 7^—, —^— n r — i °t = 


GS 
a 

G Jr) 


(15) 


(1 + iz)(l — 2v) 


a 


Then formulas (11) yield expressions, connecting param¬ 
eters of the model with characteristics of the bond: 

Bi = (1 (1 'A s —, B 2 = GSa, B 4 = ^, 
(1 + i/)(l — 2iz) a a 


B, = 


(1 - v) EJ B 2 B a 


(1 + v)(l — 2v) a 4 2 ' 


( 16 ) 

Thus in the case of short bonds, formulas (16) can be 
used for calibration of the model. 


V. RESULTS 


tially simplified. The potential energy of the bond is 
represented via vectors, rigidly connected with bonded 
particles. Forces and torques caused by the bond are 
calculated. Relations between parameters of the model 
and bond stiffnesses are established. It is shown that 
by choosing the parameters any values of longitudinal, 
shear, bending, and torsional stiffnesses of the bond can 
be fitted exactly. Therefore the model is applicable to 
bonds with arbitrary length/thickness ratio. The expres¬ 
sions connecting parameters of the model with geomet¬ 
rical and mechanical characteristics of the bond are de¬ 
rived. A method for validation of numerical implemen¬ 
tation of the model is presented. It is proposed to con¬ 
sider four types of deformations (tension, shear, bending, 
and torsion) in a two particle system and to calculate 
corresponding forces and torques. The results can be 
compared with simple analytical expressions derived in 
the present Brief Report. This approach allows to mini¬ 
mize the risk of error in numerical implementation of the 
model. 


In the present Brief Report, the model [1] for sim¬ 
ulation of elastic bonds in granular solids is substan- 
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